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ABSTRACT 

We present an algorithm for the decomposition of periodic 
financial return data into orthogonal factors of expected 
return and "systemic", "productive", and "nonproductive" 
risk. Generally, when the number of funds does not exceed 
the number of periods, the expected return of a portfolio 
is an affine function of its productive risk. 
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Preface 



This is a paper about our rtndecomp algorithm, an algorithm for decomposing 
financial return data into expected returns and principal components of risk. A 
complete listing of the algorithm appears in Appendix B. Section 5 describes exactly 
what the algorithm does. The rest of the paper is background — more or less. 

The paper is accompanied by three ancillary text files: 

• rtndecomp. m - The GNU Octave function. 

• GPLv3.txt - The GNU General Public License governing the use of the 
rtndecomp. m code. 

• AdjustedClosingPrices_2010-2011.csv - The adjusted closing prices, in tab- 
separated-value (spreadsheet) format, of 22 iShares exchange traded funds on 
the 505 market days from 2009-12-31 to 2011-12-30 inclusive. These prices are 
normalized at 100.000 on 2010-12-31. This means that the proportions in a 
notional portfolio p = (px, . . . ,^22) represent the proportions of the 22 securi- 
ties in an actual investment portfolio at the close of 2010-12-31. The security 
proportions in the same investment portfolio are typically different at the close 
of any one of the other 504 market days under consideration. 

Section 6, Examples of output, illustrates the application of the algorithm to real 
world data. All computations in this section are based on the adjusted closing prices 
in "AdjustedClosingPrices_2010-2011.csv." 
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1 The standard model 



We start with a synopsis of the "standard mean-variance portfolio selection model" 
([Markowitz(1987), pp. 3-5]) for ex post return data. 

Given an M x n matrix R = [r 1; . . . , r n ] of successive periodic returns (M returns 
for each of n securities), an investor is to choose the proportions p = [pi, . . . ,p n ] T 
invested in each security, the proportions being subject to the constraints pj > (j = 
1, . . . ,n), Y^j=\P3 = 1- We assume that the periodic returns, r p G H , of the corre- 
sponding investment portfolio satisfy the linear hypothesis 

n 

We also assume that expected periodic return is a linear function of periodic return 
or, in other words, the expected periodic return of security j is given by e 3 - = u T Vj 
for j = I, . . . , n. Here the weight vector uj G R M should satisfy Ui > (i — 1, . . . , M) 
and Y^Li^i = 1- 

Under these assumptions, the expected periodic return of the investment portfolio 
corresponding to p is 

n 

e p = Y^ e jPj = E P, ( 2 ) 
i=i 

with E — [e±, . . . , e n ] = oj t R, and the variance of portfolio return is 

n n 

v P = v ikPjPk = P T Vp, (3) 

3=1 k=l 

where the n x n covariance matrix V = [vjk] is given by 

M 

Vjk = uJiZijZik (j, k = 1, . . . n), (4) 
i=i 

the deviation or "risk" vectors Zj G R M (j — 1, . . . , n) being defined by 

z i = r i - !ivfei> (5) 
with 1m G R m representing the constant return vector of all l's. 

Caveat. If the periodic returns in R are normalized linear returns, then the normal- 
ized linear returns of each investment portfolio in the n securities satisfy the linear 
hypothesis (1) with respect to some p = [pi, . . . ,p n ] T , and all of the above statements 
follow ([Norton(2011)]). More typically, when compound periodic returns are used, 
the linear hypothesis cannot be satisfied by any p, and the arguments of this paper 
do not apply. 
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2 Geometry 



We will consider the ex post standard model from a geometric standpoint. The a 
priori weights, u?, of section 1 induce a Euclidean metric on the space of consecutive 
periodic returns, R M : 



M 



(x, y)^ = y^ y u)jXiyi, ||x|| w = ^ (x,x) w , for x,y G R M . 



(6) 



i=l 



Two return vectors x and y are orthogonal (perpendicular to each other) if (x, y) w = 0. 

The vector of all ones, 1m, is a unit vector in this Euclidean space since YlfLi u i = !• 
The expected return axis, the E-&xis, points in the lM-direction. The ^-coordinate 
of any periodic return vector r G R M , 



M 



e = (1 m ,t) u = J^. 



(7) 



is its expected return. 

Figure 1: r = z + l M e 
E 



Each periodic return vector, r, has an orthogonal de- 
composition into its (scalar) expected-return compo- 
nent, e, and its (vector) risk component, 
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U(r) 


= |z|L ^ 



z = r - l M e, 



(8) 



Z 



with expected return zero. The standard deviation of 
periodic return is simply the length or norm of the risk 
component, 

= ll z IU, (9) 
and the variance of periodic return is its square norm, 



v r 



(10) 



The covariance matrix V = [vjk] of (4) is the Gram matrix of inner products of the 
security risk vectors Z — [zi, . . . , z n ] of (5): 



(j,k = l,...,n). 



3 Linear subspaces and flats 



We are concerned with notional portfolios in n specific securities. The return vectors 
of these portfolios lie in the the linear subspace £(R) of R M spanned by the return 
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vectors, R — [r 1; . . . , r n ], of the individual securities: 

n 

L{R) = { ^rjtj : tj e R }. 

3=1 

The risk components of portfolio return vectors lie in the linear subspace £{Z) of R M 
spanned by the risk components, Z — [z^, . . . , z n ], of the r^-: 

rt 

£(Z) = {^z^:t i GR}. 

i=i 

Since the proportions of the securities in a notional portfolio must sum to 1, portfolio 
return vectors and their risk components are contained in the flats (affine subspaces 
of R M ) defined by 

n n 

W = { : I> = 1 } 

3=1 3=1 

and 

n n 

<F(Z) = 1^ = 1}, 

3=1 3=1 

respectively. We will refer to these as the R- and Z-flats. 

Finally, we will be concerned with differences in periodic return vectors, and the 
corresponding differences in their risk components, from one notional portfolio to 
another. Such difference vectors reside in the tangent spaces 

n n 

3=1 3=1 

and 

n n 

7(z) = {J2*jtr-J2 t i = °} 

3=1 3=1 

of the R- and Z-flats. 

Proposition 1. If 1m ^ 7(R), then the risk component mapping r h> z defined by 
(7) and (8) restricts to a linear isomorphism of7(R) onto 7(Z). 
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Proof. The mapping from 7(R) onto 7(Z) can be expresses as Ar i-> Az with 

n n n 

Ar = Az + 1 M ejtj, Az = z^i,-, and £j = 0. 

j=i j=i j=i 

To show that this mapping is a linear isomorphism, we need to show that Ar = Om 
whenever Az = Om- But, if Az = Om, then Ar = lMEj=i e iV And then, since 
1 M i 7(R), EJ=i ejtj = and Ar = M . ' □ 

Corollary. // 1 M ^ 7(R), then 

n n n 

j=i j=i j=i 

is a well-defined mapping of3 r (Z) onto 7(R). R is the inverse of the risk component 
mapping from 7{R) onto 7(Z). 



4 Components of portfolio risk 

The total variance of return of the periodic returns in R = [r 1 , . . . , r n ] is the sum of 
the variances of return of the individual securities: 

n n 

VT = j2vjj = j^H z jII«- ( 13 ) 

3=1 3=1 

This is a measure of the volatility of the return data as a whole, of the spread of the 
periodic returns in R = [ri, . . . , r n ] away from their expected values E — [ei, . . . , e n ]. 

Given a unit risk vector u e £(Z), the variance of return of security j in the u- 
direction is the square of the u-coordinate of its risk vector, (u,Zj)^. The total 
variance of return in the u-direction is the sum of the u-directional variances: 

n 

v u = Y(^^)l (14) 

3=1 

If U C £>(Z) is an orthonormal basis for £>(Z) (a pairwise-orthogonal set of unit 
vectors that span L(Z)), then 

\\*i\\l = Y t ( u >*i)* (J = !,••-, »)■ (15) 
ueu 
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Consequently 

n n n 

<r = E \Ml = £5>^ = E5> z ^ = ( 16 ) 

3=1 j=l u6U uSlt j=l u£tt 



For principal component analysis ([Wikipedia(2011)]) one attempts to choose the 
orthogonal basis XL (orthogonal coordinate system if you will) so that the sum of 
u-directional-total- variances on the right side of (16) decomposes or "explains" the 
total variance, v<r, in a particularly meaningful way. We are aiming for such a de- 
composition of the total variance of return in this paper. Our idea of a "particularly 
meaningful way" will be defined in this section. 



4.1 Systemic risk 

Lemma. Let z denote the point in the Z-flat that is closest to the origin: 

||z ||w = min{ ||z|| w : z G 3{Z) }. 

Then 

(z , z - z ) w = for all z G 5(Z) (17) 

Proof. Given z G ^(Z), z(t) = z + (z — z Q )t is in 3^(Z) for all t G R. By definition 
||z(t)||^ achieves its minimum value of ||z ||2, at £ = 0. Consequently 



1 d 

2 dt 



\\ z (t)\L = (zo,z - z ) w = 0. 

t=o 



□ 



Proposition 2. Let fo = ||zo|| w . Then 

Vjk = fo+Vjk with v jk = (zj - z ,z fe - z )o; (J,k = l,...,n). 



Proof. 



V jk = (Zj,Zj(.) w (11) 
= ((Zj - Zo) + Z , (Zfc - z ) + z ) w 
= (Zj - z , z fe - z ) w + (z , Z )cj 

+ (zj - z , z ) w + (z , z fe - z ) w (bilinear expansion) 

= (zj - z , z fe - z ) w + (z , z ) w + + (by the lemma) 
= Vjk + /o • 



□ 
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Corollary 1. 

n n 

vt = n/o + v T where v T = ^Vjj = \\zj - z ||£. (19) 

i=i 3=1 

This follows from (13) and (18). Then, by rewriting (18) in matrix form, we see that 
Corollary 2. 

V = fl + V where V = [v jk ] (j, k = 1, . . . , n) . (20) 

Here we adopt the convention that the sum of a scalar and a matrix is the original 
matrix with the scalar added to its every coefficient. 

We refer to fo as the systemic portfolio risk, /q is the systemic portfolio variance, 
and n/o is the total systemic variance of the system. The variance of return of any 
notional portfolio p decomposes into its systemic and nonsystemic parts: 

v P = p T Vp = f Q 2 + p T Vp. (21) 

The second equation follows from (20) and Y^=iPj = 1- 

A minimum-variance portfolio is a notional portfolio p whose variance is less than or 
equal to the variance of any other notional portfolio q with the same expected return. 
Minimum- variance portfolios play a crutial role in Markowitz's mean- variance analysis 
([Markowitz(1987)]). By (21) 

p T Vp < q T Vq if and only if p^Vp < q^Vq 

for notional portfolios p and q. Consequently, the collection of all minimum-variance 
portfolios is completely determined by the singular, nonsystemic covariance matrix 
V and the expected return matrix E. 

If fo = 1 1 Zo 1 1 cl> 7^ 0, we will take Uo = %o/fo to be the first vector in our orthonormal 
basis U for £>(Z). This is the direction of systemic risk. Every vector z G ^{Z) has 
the same Uo-coordinate, (uo, z) w = fo, as can be seen from the expansion 

(u , z) w = (u , z + (z - z Q )) w 

= (u , z ) w + (uo, z - z ) w (linear expansion) 

= (uo, u /o)w + (definition of u and (17)) 

= /o- 

We will assume until further notice that 1m ^ 7(R). Then the mapping z i— > r 
from the Z-flat onto the R-ft&t is well-defined by the corollary to Proposition 1, and 
Zo € ${Z) is the risk component of a unique ro G 3{R). We will refer to 

e = (1m, r >„ (22) 
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as the systemic portfolio return of our system. Note that eo may not be the expected 
return of any notional portfolio p, all of whose coefficients must be nonnegative. 



4.2 Productive risk 



Equation (17) shows that the tangent space 7(Z) is the orthogonal complement of z 
in &{Z). Indeed 7(Z) is spanned by the difference vectors Zj — z (j = 1, . . . , n), and 
the nonsystemic covariance matrix, V, is the Gram matrix of these difference vectors. 
We will select the remaining orthonormal basis vectors Uj (i — 1, . . . , to; m < n) from 
7(Z). Then the total nonsystemic variance of (19)) will decompose as the sum of 
the squares of the Uj-coordinates of the Zj — z , 

n m n m 

j=l i=l j=l i=l 

and V will factor as V = F T F, with the coefficients of the m x n factor matrix F 
given by 

fij = (ui,Zj-z ) u (i = 1, ...,m;j = l,...,n). (24) 



We will continue to assume that 1m ^ 7(R), so that r n- z is a bijection of onto 
7(Z), and further suppose that the n securities do not all have the same expected 
return. Under these assumptions the orthogonal projection of 1m onto 7(R) is neither 
Om nor 1m itself. 



Figure 2: The (uj., 1^) -plane 
7(E) 7{R) 













e F = 








/ 





Now set 



Ui 



Vi 



Let vi G T(i?) denote the unit vector in the di- 
rection of the orthogonal projection of 1m onto 
7(R), as shown in Figure 2. Then Vj is the di- 
rection of steepest increase of expected return in 
the i?-fiat. Changes in expected return depend 
only on changes of periodic return in the Vi di- 
rection in the sense that 



Ae =< 1 M , Ar > w 

=< 1 M , Vi > w < vi, Ar > L 



(25) 



for all Ar e 7(R). 



-m\1m, Vi)(, 



|vi - 1m (1m, vi) 



eT(Z) 



(26) 
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so that 



7T 

vx = u x cos + 1m sin <fi with < <p < — (27) 



as shown in Figure 2. 
Proposition 3. Let 



Ae = (l M ,Ar) u , 
Az = Ar - 1 M Ae, 



/or Ar G T(i2). TTien 



Ae = e F (u 1 ,Az)„ (28) 
wii/i the Ui of (26) and ep = tan as in Figure 2. 

Proof. The result follows from (25) and 
(1m,Vi) w = sin0, 

(vi, Ar) w = (ui, Ar) w cos0+ (1m, Ar) w sin0 
= (ui , Az) w cos + Ae sin 0. 

□ 



Under the assumption 1m ^ ^(-^), the change in expected return, Ae, from one 
notional portfolio to another depends only the change in risk component, Az, of 
the respective return vectors. This is a consequence of Proposition 1. Proposition 
3 now shows that such a change in expected return depends only on the change in 
the risk component in the Ux-direction. For this reason we refer to the Ux-direction 
of the Z-flat as the direction of productive risk. Changes in portfolio risk vectors in 
directions orthogonal to the Ux-direction have no effect on expected reward. Such 
changes are nonproductive in this sense. 

Corollary 1. 

e = e + e_p(ux, z) u for all r = z + l M e G 5F(i?). (29) 

Proof. Set Ar = r — r in Proposition 3. Then 

e - e = e jF (u 1 ,z - z ) w 

= ej?((ui,z) w - (ui,z ) w ) 
= e jF (u 1 ,z) tJ . 

Here (ui, zq)^ = since ui G 7(Z) and zq is orthogonal to 7(Z). □ 
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Corollary 2. 

e = e* - e i r(ui,z*) w 
/or any convenient r* = z* + l^e* G 7(R). 



(30) 



Let us now define the productive risk of the system, r 1; as 



E< u ^>s ( 31 ) 

i=l 

with rf being the productive variance. We include the middle, zq expression in this 
definition to emphasize that the productive risk is coming from the tangent space 
7(Z), which is spanned by the Zj — z . The middle expression collapses to the last 
expression because (ui,Zo) w = 0. 

4.3 Nonproductive risk 

Each notional portfolio p has a corresponding risk vector z = Zp. By (29) the 
expected return of the portfolio is completely determined by the Ui coordinate of z 
and the parameters and ep. However, the sum of the systemic and productive 
variances, /q + (u 1; z)^, is just a part of the portfolio variance. The remaining vari- 
ance is nonproductive, having no effect on the expected return of the portfolio. 

Definition 1. We now define the principal nonproductive risks, Tj > 0, and the 
corresponding principal directions of nonproductive risk, Uj G 7(Z) (||uj|| w = 1), for 
i = 2,...,m, where m is the dimension of 7(Z), the rank of V. The definition 
proceeds by induction: 

for i — 2, . . . , m, 

n 
0=1 

= max |y^(u,z 3 -)% : u G 7(Z), \\u\\ u = 1, (u k ,u) u = (A; = 1, . . . ,% - 1)| 

Remark. The = \frf are uniquely determined, and, in the generic case, when 
t~2 > f"3 > • • • > T m , the principal directions of nonproductive risk are unique upto 
multiplication by —1. We will assume this case to simplify the discussion. The 
rtndecomp algorithm presented in Appendix B makes no such assumption. 



n 
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Corollary 1. The total nonsystemic variance can be decomposed into its productive 
and nonproductive parts as 

m 

VT = t 2 1+ J2 r? (32) 

i=2 

Corollary 2. The nonsystemic covariance matrix V factors as V = F T F, where the 
coefficients of the m x n factor matrix F are given by 

fij = (ui, Zj - Zo) w = (m, Xj) u (i = 1, . . . ,m; j = 1, . . . ,n) (33) 

Corollary 3. The variance of return of any notional portfolio p can be decomposed 
into its systemic, productive, and nonproductive parts as 

n m n 

j=l i=2 j=l 

with the fij of (33). 

These three corollaries of Definition 1 follow immediately from the preceding discus- 
sion. We refer to the factor matrix F of Corollary 2 as the nonsystemic risk matrix. 



4.4 Mean- variance analysis 



Let 

denote the mean-variance mapping defined by (2) and (3). In view of the preceding 
discussion $ can be factored as 



with 



where 



$ : R n — )■ R m -> R 2 , p h-> f p h-> (e p , v p ), 



e P = e + e F x p , 

/o 2 + llf P ll 2 = /o 2 + 4 + llypl 



p = Fp. 



Xp 










F(2:m,:) 



(35) 

(36) 
(37) 

(38) 



Here F(l,:) and F(2 : m, :) denote the productive and nonproductive rows of the 
m x n nonsystemic risk matrix F, respectively. 

We are primarily interested in the image, $(A), of the notional portfolio simplex 

n 

A = { p G R" : ^Pj = 1 and pj > for j = 1, . . . , n}. 
[Markowitz(1987)] refers to this image as the obtainable EV set. 
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Remark. We apologize for the reuse of notation here. We have been using A to 
indicate a difference vector. Now A is the standard (n — l)-simplex in R n . In the 
future we hope the meaning of A will be clear by its context. 

Example 1 and the corresponding Figure 3 illustrate the factorization (35)-(37) in 
the n = 3, m = 2 case. 

Example 1. 



-4 2 
2 -2 



4 
3 





20 


-12 


-10 


F T F = 


-12 


8 


2 




-10 


2 


25 



Figure 3: Factorization of the mean- variance mapping $ 



A(l, 0, 0) 



C(0, 0, 1) 




B(0, 1, 0) 



P 

F I 

fay) 




P >-> (e,f) 

e = eo + epx 
v = fo +x 2 + ||y|| 2 



(x,y) i y {x, \y\) 



> x 




> x 



In Figure 3 the set of minimum-variance portfolios is specified by the piecewise linear 
path APQC though the simplex A. The portfolios 



and 



P = 50% A + 50% B 
Q = 60% B + 40% C 
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are called corner portfolios for obvious reasons. Equation (37) implies that the port- 
folio E of absolute minimum variance corresponds to (x,y) = (0,0). Since P and Q 
have (x, ^-representations (-1, 0) and (2.8, 0), respectively, we must have 

E=H P+ 1 Q . 
19 19 

The parameters eo, e^ > 0, and fo > are inconsequential. The set of minimum- 
variance portfolios is independent of these parameters. 

Portfolios on the piecewise linear path EQC in A are efficient: besides having min- 
imum variance for their expected return, they have maximum expected return for 
their variance. 

The x- and y-axes through the upper-left-hand simplex of Figure 3 are the preimages 
of the respective axes on the lower-left xy-plane via the mapping p i— > (x, y) = Fp. 
The image axes are perpendicular to each other, but the preimage axes are not. The 
preimage of the x-axis is the critical line of the mapping 3>|{X^=iPj = 1}- The 
derivative of &\{^2j = iPj = 1} has rank 1 along this line and rank 2 everywhere else. 
In effect $ folds the Ylj=i Pj = 1 plane over this critical line. 

Returning to the general situation let us point out that the factorization of the 
mean variance mapping in (35)-(37) leads to a natural, geometric characterization of 
minimum- variance portfolios. First note that the portfolio simplex A is mapped onto 
a convex polytope P = F(A) in (x,y)-space R m due to the linearity of p > Fp. If 
the x-axis (the Ui-axis) passes through this polytope, then every point in the inter- 
section of the x-axis and the polytope is the image of a minimum-variance portfolio 
p — simply because y p = and v p can't get any smaller than v p = /g + x p . More 
generally, let x m ; n and x max be the minimum and maximum values of the coefficients 
in the "productive" row, X = F(l, :), of F. Given x* between x m i n and x max , suppose 
y* G R m_1 satisfies 

||y*|| = min{||y|| : (x*,y) G P H {x = x*}} . (39) 

Then any portfolio p G A that F maps onto (x*,y*) (and there is at least one) is a 
minimum-variance portfolio. 

4.5 Relaxing assumptions 

Since Section 4.2 we have been assuming that the n given securities do not all have 
the same expected return; however there is no problem if the returns are identical. 
Then there is no productive risk, all nonsystemic risk is nonproductive: u 1; . . . , u m 
are the principal directions of nonproductive risk. This case is signalled by ep = 0, 
and (29) still holds with eo being the common expected return. 
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The case when 1m parallels the R-flat, when 1m G 7(R), is more problematic. This 
situation typically arises when there are more securities than periods. Then there is 
no unambiguous systemic return, eo, and no well-defined gradient of expected return, 
g = uie F G 7(Z). 

To handle the 1a/ G 7(R) case we anchor ourselves at the mean risk component, 
z = i Y^j=\ z j> w hh mean expected return, e = ^ X)j=i e i- O ur approximate gradient, 
g = X]j=i( z i — ^)9ji i s the least-squares solution of (28) in the form 

n 

e k -e = y^(z fc - z, Zj - z% (k = l,...,n). (40) 

We use this g to define the parameters 

ef=||g|| w , ui = g/e F , eo = e-(g,z) w , (41) 

for the approximate version of (29): 

e«e + ei?(ui,z) w . (42) 

Note that the definition of eo is essentially (30) with e* = e and z* = z. 

The remaining orthogonal directions of risk are defined inductively by Definition 1, 
and the corollaries of that definition continue to hold. 



4.6 Scaling output 

Up to this point expected returns, eo and E = [ei, . . . ,e n ], and risk coefficients, /o 
and F = [fy] (i = 1, . . . , m; j = 1, . . . , n), have been measured in the same percent- 
per-period units. While days or weeks may be used for computational purposes, 
annualized, percent-per-year output is usually preferred to daily or weekly percent- 
ages. 

To compensate for this preference we add a periods-per-unit-of-time parameter p 
to the periodic returns R and weights u required by our algorithm. Then, at the 
end of the computations, percent-per-period expected returns and risks are scaled to 
percent-per-unit-of-time units as follows: 

percent percent 
per-unit-of-time per-period 



E 




P 


X 


E 


e 


4- 


P 


X 


e 


e F 


4- 


Vp 


X 


e F 


F 


4- 


Vp 


X 


F 


fo 


4- 


Vp 


X 


fo 
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The idea behind this scaling is statistical. Assume, for example, that daily returns 
are independent random variables from one market-day to the next and there are 
(typically) p = 252 market-days per year. The annual return is the sum of p daily 
returns; so the expected value of annual return is p times the daily expected value. 
This accounts for the p multipliers above. The variance of annual return is p times the 
daily variance due to the independence assumption, but risk or standard deviation is 
the square root of variance; consequently yjp~ is the appropriate multiplier of F and 
/ . Finally ep is the rate of change of expected return to risk; so yfp = p/sfp is the 
appropriate multiplier. 



5 The rtndecomp function - arguments and relation- 
ships 



The GNU Octave listing of the rtndecomp function appears in Appendix B. In this 
section we give the function header and describe its arguments. The relationships 
between the output arguments were derived in the last section. 

function: [E, F, f , e , ep] = rtndecomp (R, a;, p) 

purpose To decompose financial return data into orthogonal risk-factors. 



input 

R 



output 

E 
F 

fo 

e F 



- M x n matrix of periodic returns. 

- M-vector of positive weights or a scalar. 

If u? is a scalar or if R is the only input argument, 
then u) defaults to 

Ui = 1/M for i = 1,...,M. 

- periods per unit time, (p > 1, default: p = 1) 
e.g., p = 252 market-days per year. 

- 1 x n matrix of expected returns. 

- m x n matrix of risk coefficients. 
rank(F) = m unless F = zeros(l,n). 

- systemic risk, (fo > 0) 

- systemic expected return. 

- expected return per unit of productive risk. 
{&f > 0; if Cf = there is no productive risk) 



global output 
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eflag - true if and only if a nonzero, constant M-vector is 

parallel to the returns flat, J(R), or, said another way, 
if and only if 1a/ G 7(R). 

relationships 

1) E = pu> T R. 

2) The n x n covariance of returns matrix, V, is given by 

V = Z T diag(up) Z, where Z = R- l M u T R. 

3) V = /o + F T F. (the scalar /q is added to each coefficient of F T F) 

4) E = Cq + epX unless eflag is true. 

Here X = F(l,l : n) denotes the first row of F, and, in the equation, Cq 
is added to each coefficient of epX. When eflag is true, the equation is an 
approximation, but 

mean(E) = e + e^mean(X) 
remains true. 

5) v T = J2U V(j.j) = up, + rl + £r= 2 7f , with n = y/F(i, 1 : n)F(i, 1 : n) T 
for i = 1, . . . ,m. The right-hand side of the equation for vt represents the 
decomposition of total variance into systemic, productive, and nonproduc- 
tive parts (though r\ is nonproductive if ep = 0). 

6) T2 > . . . > T m > 0. (principal nonproductive risks) 

7) F(2 : m, 1 : n)F(2 :m,l:n) T = diag(r|, . . . , r^). 

8) If e F = and F ^ [0, . . . , 0], then n > . . . > r m > 
and F(l : m, 1 : ra)F(l : m, 1 : n) T = diag(r^, . . . , r^J. 



6 Examples of output 
6.1 Five large ETFs 2010 

In our first example we apply the rtndecomp algorithm (listed in Appendix B) to 
2010 daily returns from the five iShares exchange traded funds (ETFs) 

1. IEF - iShares Barclays 7-10 Year Treasury Bond Fund 

2. IWB - iShares Russell 1000 Index Fund 

3. IWM - iShares Russell 2000 Index Fund 

4. EFA - iShares MSCI EAFE Index Fund 

5. EEM - iShares MSCI Emerging Markets Index Fund 

Figure 4 shows the growth of these securities in 2010. These are plots of adjusted 
closing prices against time. The adjusted closing prices are normalized at 100 on 2010- 
12-31; thus notional portfolios specify the closing proportions of actual investment 
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portfolios at the end of 2010 ([Norton(2011)]). If ex post analysis were to deem a 
certain notional portfolio p as optimal, one would buy 



100 Pj/cij 



(43) 



shares of fund j, per $100 invested, to invest in the optimal portfolio at the end of 
2010. Here aj is the 2010-12-31 closing price of fund j (j = 1, . . . , n; n = 5). 

Figure 4: 2010 adjusted closing prices of five large ETFs 
prices normalized at 100 on 2010-12-31 

110 



100 




90 



80 



70 



2009-12-31 



6.1.1 Algorithm input 



The graphs in Figure 4 correspond to a 253 x 5 matrix A of adjusted closing prices 
for these 5 ETFs on the 253 market days from 2009-12-31 thru 2010-12-31 inclusive. 
The normalized, linear, daily returns for these funds are given by the 252 x 5 matrix 
R = AA, where A is the 252 x 253 difference operator 



A 



1 



-1 











-1 1 



We won't use all of R in this example, just the last M = 200 rows, just the last 200 
market-day returns — which correspond to the colored-background portion of Figure 
4. 
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In addition to the M x n return matrix i?, the rtndecomp algorithm requires an 
M- vector of weights, lo, and a scaling factor, p, that specifies the number of periods 
per unit time. We use p = 252 market-days per year throughout this paper. 

To see how the weights affect output we will consider two systems of weights: gray, 
uniform weights, where each of the M = 200 market days has the same importance, 
bJi = 1/200, and the more colorful late-heavy weight system pictured in Figure 5. 



The colors of the late-heavy system corre- 
spond to the colored regions of Figure 4. 
The returns of the last 30 market days of 
2010 are weighted by 1/140 each; returns 
for the first 70 days count half as much as 
these, or 1/280 per day; and the weights for 
the middle 100 days, 



Figure 5: Return weights for 200 
market days. Late-heavy weights 
colored. 



to 



z/101 



280 



for 



100, 



1/140 

1/200 
1/280 



increase uniformly between these two ex- 
tremes. The sum of these late-heavy weights 
is the sum of the yellow, green, and blue ar- 
eas in Figure 5: 



uniform 



> i 



70 



170 200 



1 1 

70 x + 100 x - x 

280 2 



1 1 \ 1 

+ + 30 x 

280 140 J 140 



as required. 



The idea behind late-heavy weighting is simple. Think of investing at the end of 
a 200 market-day period. The recent performance of a group of securities may be 
more important than their performance further back in the past — as a predictor of 
their near-future performance. Thus an investment analysis might weight the recent 
performance more heavily. 



6.1.2 Algorithm output 

The output arguments of the rtndecomp function 

[E, F, f , e , e F ] = rtndecomp( J R, u, p) (44) 

were described in the Section 5. Tables 1 and 2 show the output when the function 
is applied to the 200 x 5 matrix of daily returns described above. The late-heavy 
weights of Figure 5 were used for Table 1 and uniform weights for Table 2. 



17 



Table 1: Decomposition of return data - 5 large ETF universe 
last 200 market-days of 2010 - late-heavy weights 



fund 


IEF 


IWB 


IWM 


EFA 


EEM 




E 


2.86 


18.77 


27.12 


13.84 


21.79 


V T 




-2.79 


7.80 


13.36 


4.52 


9.81 


364 


24.8% 


F 


4.14 


-11.81 


-13.48 


-19.57 


-16.81 


1004 


68.4% 




2.30 


-2.15 


-5.06 


0.01 


6.12 


73 


5.0% 




-4.14 


2.02 


-2.26 


-0.89 


0.40 


27 


1.9% 




47 


209 


391 


404 


416 


1468 


100% 


3.2% 


14.2% 


26.6% 


27.5% 


28.4% 


100% 






with /o 


= 5.07, 


e = 7.05, 


and ep 


= 1.502 







Table 2: Decomposition of return data - 5 large ETF universe 
last 200 market-days of 2010 - uniform weights 



fund 


IEF 


IWB 


IWM 


EFA 


EEM 




E 


7.79 


11.62 


17.47 


8.12 


17.27 


V T 




0.00 


3.02 


7.62 


0.27 


7.46 


123 


7.4% 


F 


5.41 


-15.08 


-18.63 


-20.92 


-19.15 


1408 


85.3% 


-1.74 


2.31 


6.80 


-2.94 


-5.71 


96 


5.8% 




-3.33 


2.09 


-1.48 


-2.17 


1.22 


24 


1.4% 




43 


246 


453 


451 


456 


1651 


100% 


2.6% 


14.9% 


27.5% 


27.3% 


27.7% 


100% 






with /o 


= 4.86, 


e = 7.78, 


and ep 


= 1.271 







The E, eo, F, and /o coefficients in these tables are in percent-per-year units, the 
slope ep is unitless, and the nonsystemic variance totals in the Vt sections are in 
percent-per-year-squared units. 

As described in Section 4, the coefficients of each risk matrix F are the coordinates 
of the nonsystemic risk components of the individual securities with respect to an 
orthonormal basis, {iii,...,u m }, for the tangent space of the Z-fiat. The sum of 
the squares of these coefficients is the total nonsystemic variance of the system. The 
Vr-row of either table shows how this nonsystemic variance is distributed among the 
individual funds. As one might expect, the nonsystemic variance of the bond fund, 
IEF, is substantially less than that of any of the equity funds — under either system 
of weights. 

The green row, X = [xi, . . . ,x n ] = F(l, :), of each F contains the productive risk 
coefficients, the Ui-coordinates, of the security risk vectors. Changes in productive 
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risk from one notional portfolio to another produce corresponding changes in expected 
return. If there is no change in productive risk, there is no change in expected return. 



The orange rows, F(2 : m, :), of each F contain the risk coefficients in the principal 
directions of nonproductive risk. These rows are pairwise orthogonal (up to roundoff 
error). The first orange row, Y — [yi, . ■ ■ ,y n ] — F(2, :), represents the most significant 
or major direction of nonproductive risk. It contains the u 2 -coordinates of the security 
risk vectors. The Vr-column of each table shows how the total nonsystemic variance 
is decomposed into productive and nonproductive components. 

In addition to its nonsystemic variance each fund has a systemic variance of /q so 
that the total variance of the system is decomposed into its systemic, productive, and 
nonproductive parts as 

VT = nf + EE 4" (45) 

j=l i=2 j=l 

Table 3 shows the total variance decompositions corresponding to Tables 1 and 2. 



Table 3: Decomposition of total variance 
5 large ETFs - last 200 market-days of 2010 





late-heavy 
weights 


uniform 
weights 


systemic variance (ra/g) 
productive variance (Yl x2 ) 
major nonproductive variance y 2 ) 
other nonproductive variance 


128 8.0% 
364 22.8% 
1004 62.9% 
100 6.3% 


118 6.7% 
123 7.0% 
1408 79.6% 
120 6.8% 


total variance 


1596 100.0% 


1769 100.0% 



6.1.3 The XE-plane 

We will let 

n 

a; = Xp = F(l,:)p = Ete (46) 

i=i 

denote the productive risk coordinate of z = Y^j=i z jPj anc ^ 

n 

e = Ep = J2e jPj (47) 

be the corresponding expected return coordinate. Figure 6 shows the graphs of 

e = eo + epx, (48) 
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corresponding to the late-heavy and uniform weight systems. The plotted security 
points realize the X and E rows of Tables 1 and 2. The grid scale is 5 x 5. 



Figure 6: expected return as an 
affine function of productive risk 
last 200 market days of 2010 




late-heavy weights uniform weights 



These are 2-dimensional slices of the 
6-dimensional ZE-sp&ces, £(Z) x R, 
corresponding to the two weight systems 
we are considering. Each periodic return 
vector, r = z + l^e, corresponds to a point 
(z, e) in ZE-sp&ce. Figure 6 shows the 
orthogonal projection, (x,e), x = (ui,z) w , 
of these points onto the respective XE- 
planes. 

In each picture the blue segment connect- 
ing the security points is the projection of 
the portfolio polytope in ZE-sp&ce onto 
the XE-pl&ne. The line through this blue 
segment is the projection of the entire R- 
flat. 



Figure 6 shows how the relationship between productive risk and expected return can 
vary with weight system. Productive risk accounts for 22.8% of the total variance 
under the late-heavy system, but only 7.0% under the uniform system. This difference 
shows up in the extra width of the late-heavy picture. 

Looking at the vertical spread of Figure 6 one notices that the expected returns of the 
stock funds are higher and the bond fund return lower under the late-heavy system. 
This is because the negative returns of the stock funds in the first 70 days of the 200 
market-day sample count more in the uniform system, and the negative returns of 
the bond fund in the last 30 days count more in the late-heavy system. (Figure 4) 



6.1.4 The eflag flag 

When the matrix equation E = + epX holds exactly, the global variable eflag is 
or false. This is typically the case when there are many more periods than funds, 
as in the five-fund, 200-market-day examples just considered. Figure 6 displays this 
relationship graphically. 

Rather than increase the number of funds to illustrate the eflag = true condition let 
us decrease the number of periods from M = 200 market days to M = 3 quarters. 
Now p = 4 (quarters per year), and we will use late-heavy weights, u> = [ou\, U2, 0O3] T , 
comparable to those in Figure 5. 

Here is the complete data. 
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Table 4: Quarter-ending adjusted closing prices 



date 


IEF 


IWB 


IWM 


EFA 


EEM 


2010-03-31 


92.925 


91.157 


85.766 


93.619 


87.102 


2010-06-30 


100.196 


80.605 


77.316 


79.141 


77.686 


2010-09-30 


104.498 


89.958 


85.884 


93.452 


93.194 


2010-12-31 


100.000 


100.000 


100.000 


100.000 


100.000 



Table 5: Quarterly returns and late-heavy weights 



quarter 


IEF 


IWB 


IWM 


EFA 


EEM 


weights 


2 


7.271 


-10.552 


-8.450 


-14.478 


-9.416 


2/9 


3 


4.302 


9.353 


8.568 


14.311 


15.508 


3/9 


4 


-4.498 


10.042 


14.116 


6.548 


6.806 


4/9 



quarterly returns = adjusted closing price differences 



Consider the return vector matrix, R = [r 1; . . . , r 5 ], in Table 5. It is easy to see that 
the return-flat tangent space 

I 3=1 3=1 ) 

is all of R 3 ; in particular, the constant return vector, 1 3 , is contained in 7(R). This 
is the eflag = true condition implying that equation (48) is not exact. 

The rtndecomp output corresponding to the data of Table 5 is shown in Table 6. 
The yellow row shows the approximate e-values that result from applying equation 
(48) to the x-values of the green, productive risk row. 



Table 6: Decomposition of return data - 5 large ETF universe 
last three quarters of 2010 - late-heavy weights 



fund 


IEF 


IWB 


IWM 


EFA 


EEM 




E 


4.20 


20.94 


29.01 


17.85 


24.41 




approx 


4.84 


24.15 


26.64 


21.75 


19.03 


F 


-10.13 


12.42 


15.34 


9.63 


6.45 


626 


42.2% 


-0.78 


-11.44 


-8.37 


-19.06 


-17.07 


856 


57.8% 


V T 


103 


285 


305 


456 


333 


1483 


100% 


3.2% 


14.2% 


26.6% 


27.5% 


28.4% 


100% 






with f = 


0.005, e 


= 13.51, 


and ep 


= 0.856 
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Figure 7 shows the projection, (x,y,e) h> (x, e), 
of the Table 6 data onto the XE-p\axie. Here y is 
the nonproductive risk variable represented by the 
orange row of F. The blue polygon is the image of 
the portfolio polyhedron. 

The line through the blue polygon is the graph 
of the (approximate) expected return function 
e = eo + epx (48). The rtndecomp algorithm 
guarantees that this graph passes through the mean 
XE-point, (x,e) = (6.74, 19.28) in this example. 



Figure 7 is comparable to the late-heavy side of Fig- 
ure 6. The e- values have the same order in both 
pictures, but the x-values of the middle funds, EFA, 
IWB, and EEM, are permuted from one picture to the 
other. 

6.1.5 The IF- and EV- planes 

Let us now return to the late-heavy weight output in Table 1. Given a portfolio p let 

x = x p = Xp = F(l,:)p, 
y = y p =Yp = F(2,:)p, 

denote the productive and major nonproductive risk coordinates of p (or really of 
z p = Zp), respectively, and let 

y = y P = F(2 : m, :)p, m = 4, 

denote the full vector of nonproductive risk corresponding to p. As noted in Section 
4.4, 

e = e p = e + e F x p , (36) 

v = v p = /o + xp + ||y p f, (37) 

Figure 8 shows the images, (X, F)(A) and $(A), of the portfolio simplex A in the 
XY- and £V-planes, respectively. The x and y grid lines are 5 units apart. Since 
the (x,y)-tuples are coordinate vectors with respect to an orthonormal basis, the 
XF-image is the perpendicular projection of the four-dimensional polytope, F(A.), 
onto the XY-plane. 

In Figure 8 the ^/-coordinates of the securities are actually the negatives of the Y 
coefficients in Table 1. This sign change makes the comparison of the XY and EV 



Figure 7: expected return 
as an approximate function 

of productive risk 
last three quarters of 2010 

e 



IWM 




— ' — 1 — | — 1 — 1 > x 

late-heavy weights 
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Figure 8: Planar representations of return data - 5 large ETFs 
last 200 market-days of 2010 - late-heavy weights 



XY-p\ane 
grid = 5x5 



EV-p\ane 

e = e + epx 



y 



v = fo +x 2 + ||y 



EFA 



EEM 




EFA 



IEF 



-0.46 7.05 14.57 22.08 29.59 -0.46 7.05 14.57 22.08 29.59 
p E = 69.0% IEF + 31.0% IWB 

images more natural, but it has no effect on our analysis — a principal direction of 
nonproductive risk is, at most, determined up to a reflection through the origin. 
On the other hand, it is important to note that the stock funds are all in the first 
quadrant and the bond fund, IEF, is in the third quadrant of XY-side of Figure 8. 
This corresponds to the fact that the stock funds are positively correlated with each 
other and negatively correlated with the bond fund. This is also why the bond fund, 
IEF, is a component of every minimum-variance portfolio other than single security 
portfolio of maximum expected return, IWM. 

The solid black path through either image corresponds to the set of minimum- variance 
portfolios. As noted at the end of Section 4.4, a minimum-variance portfolio at a 
particular x = x* must minimize the value ||y|| on the polytope F(A.) fl {x = x*}. 
Points on the dotted path in either image approximate this criterion. They correspond 
to portfolios that minimize \y\ rather that ||y||. 

The point E in either figure is the image of the portfolio pe of absolute minimum 
variance, pe is an efficient portfolio. The solid black path to the right of E is the 
image of the other efficient portfolios. These efficient portfolios, in total, make up 
the piecewise linear path in A that goes from pe through the "corner portfolios" of 
Table 7 to IWM. The corner portfolios show up as the corners above the x-axis in the 
XF-image of the minimum-variance path. 
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Table 7: Optimal portfolio paths - 5 large ETF universe 
last 200 days of 2010 - late-heavy weights 





minimum- 1| y || 


(efficient) path 


minimum- y 


VjCXi Lll 




Pe 


corner 
portfolios 


IWM 


XE 

portf 


corner 
portf 


IWM 


IEF 


U.UjU 


0.618 


0.550 


n 


777 

U. i it 


U. ( \JO 


n 


IWB 


0.310 


0.382 

















IWM 








0.450 


1.000 


0.187 


0.235 


1.000 


EFA 














0.036 








EEM 























X 


0.49 


1.26 


4.47 


13.36 


0.49 


1.00 


13.36 


e 


7.79 


8.94 


13.78 


27.12 


7.79 


8.56 


27.12 


a 


5.69 


5.88 


8.48 


20.41 


6.33 


6.38 


20.41 


avg e 




17.45 






17.45 




rms a 




12.70 






12.74 





Table 7 also shows the xe- and corner porfolios of the minimum- \y\ path over the 
efficient x-range from Xe = 0.49 to x max = 13.36. The minimum- \y\ path and the 
efficient path are exactly the same from x = 4.47 to x max , but the paths differ between 
Xe and x = 4.47, the most substantial cr-differences occurring near Ie- 

The average value of e over the two portfolio paths in Table 7 is just the average of 
the end values, 7.79 and 27.12. On the other hand, the average variance, 

avgt> = / v(x)dx, 

^max •'C'Ei J xe 

and the root-mean-square risk, rms a = y/avgv, depend on the whole path. 

Remark. Throughout this paper we use Markowitz's Critical Line Algorithm as 
described in [Niedermayer and Niedermayer(2006)] to compute minimum-variance 
paths through portfolio simplices. 

6.2 Eighteen emerging markets ETFs 2010 

Now let us consider a larger universe of securities — the 18 iShares emerging markets 
ETFs that existed throughout 2010 

1. BKF - iShares MSCI BRIC Index Fund 

2. ECH - iShares MSCI Chile Investable Market Index Fund 

3. EEM - iShares MSCI Emerging Markets Index Fund 

4. EMIF - iShares S&P Emerging Markets Infrastructure Index Fund 
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5. 


EPU 


iShares 


MSCI All Peru Capped Index Fund 


6. 


ESR 


iShares 


MSCI Emerging Markets Eastern Europe Index Fund 


7. 


EWM 


iShares 


MSCI Malaysia Index Fund 


8. 


EWT 


iShares 


MSCI Taiwan Index Fund 


9. 


EWW 


iShares 


MSCI Mexico Investable Market Index Fund 


10. 


EWY 


iShares 


MSCI South Korea Index Fund 


11. 


EWZ 


iShares 


MSCI Brazil Index Fund 


12. 


EZA 


iShares 


MSCI South Africa Index Fund 


13. 


FCHI - 


iShares 


FTSE China (HK Listed) Index Fund 


14. 


FXI 


iShares 


FTSE China 25 Index Fund 


15. 


ILF 


iShares 


S&P Latin America 40 Index Fund 


16. 


INDY - 


iShares 


S&P India Nifty 50 Index Fund 


17. 


THD 


iShares 


MSCI Thailand Investable Market Index Fund 


18. 


TUR 


iShares 


MSCI Turkey Investable Market Index Fund 



Our input to rtndecomp will be the 200 x 18 matrix R = [ry] of normalized, linear, 
daily returns for the 18 emerging markets funds listed above, over the last 200 market 
days of 2010. The returns are normalized on 2010-12-31 — they are daily adjusted- 
closing-price differences divided by 2010-12-31 adjusted closing prices. We will stick 
to the late-heavy weights us of Figure 5 and use p = 252 market-days per year. 

The 17 x 18 risk matrix F corresponding to this example is not displayed, but Table 8 
summarizes how the total variance of return is decomposed by fo and F, and Figure 
9 shows the XY and EV planar representations of the rtndecomp output. 

Table 8: Decomposition of total variance - 18 emerging market ETFs 
last 200 market-days of 2010 - late-heavy weights 



systemic variance (riffi) 
productive variance (]TJ x 2 ) 
major nonproductive variance i^2,y 2 ) 
other nonproductive variance 


1131 13.9% 
1090 13.4% 
4427 54.3% 
1512 18.5% 


total variance 


8160 100.0% 



As in Section 6.1.5, E is the image of the efficient portfolio, pe, of absolute minimum 
variance. All efficient portfolios from this 18 ETF universe are made up of the four 
funds circled in red. These four funds have the least risk (a) of the eighteeen, and 
their expected returns are among the highest. This is an unusual situation — where 
risk and return seem to be inversely related. 

As in Section 6.1.5, the solid black path in either picture corresponds to the set of 
minimum- variance portfolios. The path of minimum- \y\ portfolios is dashed. The 
faint interior lines are two-security-portfolio paths. 
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Figure 9: Planar representations of return data - 18 emerging markets ETFs 
last 200 market-days of 2010 - late-heavy weights 



XF-plane 
grid = 5x5 



.EV-plane 

e = eo + epx 



y 



v = fZ + x* + \\yf 



y = 25 



A 




TUR 



EPU 



■> e 



-2.71 15.57 33.86 52.14 -2.71 15.57 33.86 52.14 
p E = 39.5% ECH + 13.0% EPU + 41.6% EWM + 5.9% EWT 

The efficient, minimum- ||y|| portfolio at x — x ECH = —3.33, 

p = 41.5% ECH + 24.0% EPU + 29.3% EWM + 5.1% EWT, 

has expected return e = 39.97 and risk a = 12.71. On the other hand, the single 
security ECH is the minimum- \y\ portfolio at this value of x (and e), but the risk of 
ECH is (T E ch = 14.46, or 13.79% more than the efficient value. Apparently, in this case, 
the minimum-||/| portfolio is not a good approximation of the minimum- ||y|| portfolio. 
This is apparent in Figure 9. 

6.2.1 The efficient four 

Let us restrict our attention to the four emerging market funds, ECH, EPU, EPU, EPU, 
that make up the efficient portfolios of Figure 9. Figure 10 shows how these funds 
grew in 2010, and Table 9 shows the output of rtndecomp restricted to their returns 
over the last 200 market days of 2010. 
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Table 9: Decomposition of return data - 4 emerging market ETFs 
last 200 market-days of 2010 - late-heavy weights 



fund 


ECH 


EPU 


EWM 


EWT 




E 


39.97 


52.07 


30.71 


36.20 


V T 




1.45 


8.51 


-3.94 


-0.74 


91 


31.0% 


F 


-5.98 


5.01 


3.04 


7.55 


127 


43.6% 




-3.44 


3.58 


3.04 


-6.33 


74 


25.4% 


V T 


50 


110 


34 


98 


292 


100% 


17.0% 


37.8% 


11.7% 


33.5% 


100% 





with f = 12.62, e = 37.47, and e F = 1.717 



The XY and EV planar representations of Table 9 are shown in Figure 11. 

The minimum-||y || portfolio path in Table 10 is efficient in either the current four- 
fund universe or in the original eighteen-emerging-market-fund universe. On the other 
hand the minimum |y|-path is based on the four-fund XY representation in Figure 
11. This minimum- \y | path is an extremely close approximation of the efficient path, 
with the maximum cr-difference of less than 0.2% occurring at x = Xe = 0. 
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Figure 11: Planar representations of return data - 4 emerging markets ETFs 
last 200 market-days of 2010 - late-heavy weights 
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Table 10: Optimal portfolio paths - 4 emerging market ETF universe 
last 200 days of 2010 - late-heavy weights 





minimum- y 


(efficient) path 


minimum- y 


path 




Pe 


corner 
portfolios 


EPU 


x E 
portf 


corner 
portf 


EPU 


ECH 


0.395 


0.464 


0.455 





0.371 


0.456 





EPU 


0.130 


0.503 


0.545 


1.000 


0.156 


0.544 


1.000 


EWM 


0.416 











0.473 








EWT 


0.059 


0.033 

















X 


0.00 


4.93 


5.30 


8.51 


0.00 


5.29 


8.51 


e 


37.47 


45.94 


46.57 


52.07 


37.47 


46.55 


52.07 


a 


12.62 


13.55 


13.70 


16.42 


12.64 


13.69 


16.42 


avg e 




44.77 






44.77 




rms a 




13.73- 






13.73+ 
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6.3 Relative risk decomposition 



The risk, er, of an individual fund or portfolio depends only on its periodic returns. 
However, the systemic, productive, and nonproductive components of this risk depend 
on the universe of funds in which the fund or portfolio resides. Table 11 illustrates 
this dependence with a fund and a portfolio from the universes we have considered. 
In this table each total risk, a, is the square root of the sum of the squares of its four 
component risks. 



Table 11: Decomposition of risk relative to fund universe 
last 200 days of 2010 - late-heavy weights 





decomposition 


decomposition 




of EEM 


of 


Pe 




18 fund 


5 fund 


18 fund 


4 fund 




universe 


universe 


universe 


universe 


expected return (e) 


21.79 


21.79 


37.47 


37.47 


systemic risk (/o) 


7.93 


5.07 


7.93 


12.62 


productive risk ( \x ) 


8.30 


9.81 


4.01 


0.00 


major nonproductive risk ( \y |) 


17.37 


16.81 


7.83 


0.00 


other nonproductive risk 


2.96 


6.13 


4.37 


0.00 


total risk {a) 


21.03 


21.03 


12.62 


12.62 



6.4 2011 results 

So far we have restricted our examples to the last 200 market days of 2010 with 
adjusted closing prices and returns normalized at the closing prices of that year. 
Late-heavy weights have been emphasized with the idea that a strong performance 
in the latter part of 2010 should carry over into 2011. 

This did not turn out to be the case. Figures 12 and 13 show how the five large ETFs 
of Section 6.1 and the four emerging market ETFs of Section 6.2.1 performed over 
2011. These are graphs of daily adjusted closing prices. Again the prices have been 
normalized at 100 on 2010-12-31 so that notional portfolio proportions correspond to 
investment portfolio proportions at 2010 closing prices. 
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Figure 12: 2011 adjusted closing prices of five large ETFs 
prices normalized at 100 on 2010-12-31 
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Figure 13: 2011 adjusted closing prices of four emerging market ETFs 
prices normalized at 100 on 2010-12-31 
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Tables 12 and 13 show the output of rtndecomp as applied to this data. Here we 
have used the full 252 markets-days of returns with uniform weighting. Now the 
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expected returns are the total returns of the respective securities over the whole of 
2011. 



Table 12: Decomposition of return data - 5 large ETF universe 
the 252 market-days of 2011 - uniform weights 



fund 


IEF 


IWB 


IWM 


EFA 


EEM 




E 


15.65 


1.24 


-4.43 


-12.23 


-18.79 


V T 




1.59 


-6.29 


-9.39 


-13.66 


-17.26 


615 


20.8% 


F 


-6.91 


21.05 


27.09 


22.92 


21.30 


2203 


74.6% 


-1.80 


0.21 


-5.70 


7.11 


-1.19 


88 


3.0% 




-1.84 


1.27 


-3.35 


-2.36 


4.96 


46 


1.6% 


Vt 


57 


484 


866 


768 


777 


2952 


100% 


1.9% 


16.4% 


29.3% 


26.0% 


26.3% 


100% 






with f 


= 4.62, e 


= 12.74, and e F 


= 1.827 







Table 13: Decomposition of return data - 4 emerging market ETF universe 
the 252 market-days of 2011 - uniform weights 



fund 


ECH 


EPU 


EWM 


EWT 




E 


-26.45 


-21.78 


-2.60 


-21.78 


V T 




-10.85 


-8.00 


3.70 


-8.01 


260 


28.8% 


F 


4.02 


-19.66 


3.58 


7.08 


466 


51.7% 




9.09 


-1.70 


-0.84 


-9.46 


176 


19.5% 


V T 


217 


453 


27 


204 


901 


100% 


24.0% 


50.3% 


3.0% 


22.6% 


100% 





with f = 19.88, e = -8.66, and e F = 1.639 



It is quite easy to show that the expected expected returns in Tables 12 and 13 are, 
in fact, total returns for the year. In the general the expected return of a security is 
given by 

M 

e = p^UiTi (49) 
i=i 

with the Ti (i = 1, . . . , M) being the successive periodic returns. When p = M, 
Ui = 1/M, and = (cij — di_i)/a , with the aj (i = 0, 1, . . . ,M) being successive 
adjusted closing prices for the security, equation (49) simplifies to 

e = ^ - 1, (50) 

which is, essentially by the definition of adjusted closing prices ([Norton(2010)]), the 
total return of the security over the M periods. 
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7 Summary 



We have described an orthogonal decomposition of the space of ex post periodic re- 
turns for a given universe of securities. The risk space, which is orthogonal to the 
expected return axis, is decomposed into systemic, productive, and principal nonpro- 
ductive dimensions. Our rtndecomp function accomplishes this decomposition. A 
technical discussion and listing of this algorithm is given in the appendix. 

In the examples of Section 6 we have emphasized the two-dimensional XF-projection 
of periodic return data. The minimum- \y \ path through the portfolio simplex is easily 
obtained from the XV-projection of the data. The minimum- \y \ path can very closely 
approximate the minimum-variance path for a small universe of securities. 

In the future we hope to develop a minimum-variance algorithm of the form 

function: P = minvar(i?, F) 
with [E, F] = rtndecomp (R, lj, p). 

Here P = [pi, . . . , p np ] would contain successive corner portfolios of the minimum- 
variance path through A and include single security portfolios at either end. 

Given [E, F] = rtndecomp (R, u>, p) and such a minvar function, the minimum-||/| 
path would be given by P 2 = minvar(£', F(l : 2,:)). More generally, the piecewise- 
linear paths through A determined by the portfolios in P^ = minvar^, F(l : k, : 
)) (k = 2, ...,m) would approximate the minimum- variance path P = P m with 
successively better approximations; moreover each portfolio in Pk would contain at 
most k securities. 
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A The flow of the rtndecomp algorithm 



Section B gives the complete GNU Octave listing of the rtndecomp algorithm. In 
this section we describe the ideas behind specific sections of the listing. Certain 
simplifying assumptions have been made in the interest of clarity. For example we 
only consider p = 1 periods per unit time. Scaling for different p, e.g. p = 252 
market-days per year, occurs at the end of the algorithm, as described in Section 4.6. 

This section is arranged in blocks. Each block summarizes a specific section of the 
rtndecomp code. 



Initial setup (rtndecomp: 105 - 111) 
E = uj t R- # uj > M , lj t 1 m = w = 1. 

/3 = v^; # ||/3|| 2 = /3 T /3 = 1. 

Z = diag(/3)i? — f3E; (3 T Z = 0, V = Z T Z (covariance matrix). 

The linear isometry 

r i — y x — diag(/3)r, f3 = 

converts the w-metric into the standard, sum-of-squares-metric for R A/ : (r, r) w = 
x T x. The risk vectors Zj in the M x n matrix Z = [z 1; ...,z n ] are the isometric 
images of the risk vectors Zj of (5); now the covariance matrix V = [vjk\ of (4) is 
given by V = [zjz k ] = Z T Z. 



QR factorization of Z (rtndecomp: 113 - 121) 
Z = QF; # compact QR factorization, 

(now F ~ [zi,z 2 , . . . ,z n ]) 

The Octave code for the QR factorization actually reads 

[Q, F, J] = qr(Z, 0); 
The output consists of an M x n matrix Q with orthonormal columns, an upper- 
triangular, n x n matrix F, and a permutation, J, of the index sequence [1, . . . ,n]. 
These matrices satisfy Z(:, J) = QF. 

This is QR factorization with column pivoting. At the end of the actual rtndecomp 
algorithm, F-columns are returned to the initial order with the replacement F : = 
F(:, J^ 1 ). In this description we assume that M > n and that Z is not rank deficient. 
Thus we can skip column pivoting and start with the QR factorization Z = QF. 

The columns of Q make up an orthonormal basis for the range of Z in R M . The 
columns fi,...,f n of the upper-triangular F are the coordinate vectors of the 
Zi,...,z n in Z with respect to this basis. We signify this situation with the no- 
tation F ~ [zi, . . . , z n ]. It simply says that "F is the matrix of coordinate vectors for 



33 



[zi, . . . , z n ] with respect to some orthonormal basis for the range of [z 1; . . . , z„]." Then, 
regardless of the orthonormal basis Q, zjz^ = $JQ T Q ffc = ffffc for j, k = 1, . . . , n, 
since Q T Q = I n . 



Z-fiat tangent space (rtndecomp: 158 - 163) 
F(l,j):=F{l,j)-F{l,l); (j = 2,...,n) 
(now F ~ [zi, z 2 - zi, . . . , z n - zi]) 

The matrix F is still upper triangular, but now the z-vectors corresponding to the 
second through the last columns of F span the Z-flat tangent space, 7(Z). 



Hessenberg QR via Givens (rtndecomp: 165 - 175) 
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±/o 

















(A.l) 



(Still F rsj [ Zl , z 2 - Zi, . . . , Z n - Zi]) 

Here we apply a sequence of n — 1 Givens rotations, G\, . . . , G n _i, to zero the subdi- 
agonal elements of the upper Hessenburg submatrix F(:, 2 : n). Then 

F := Gl_ x ■ ■ ■ GjF and Q:=QGf- 

with QF = [zi, z 2 — Zi, . . . , z n — z\] at either end of the sequence. This process is 
described in Section 5.2.4 of [Golub and Van Loan(1989)]. 

Looking at the last row of the final F in (A.l) we see that q n is orthogonal to the 
Z-flat tangent space: q^(zj — z x ) = for j = 2, . . . ,n. It follows that z = q n /ni 
is the point on the Z-flat that is closest to origin, with f Q = |/ nl | = ||z || being the 
systemic risk of the system. 



Extract systemic risk f (rtndecomp: 177 - 189) 
fo — I /nil; 

discard the last row of F: 
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m 



x 
x 
x 
x 

L ±/o 

n — 1; 



(now F 



xxx 
xxx 
x x — • 
x 


the number of rows of F) 

Zi, . . . , Z n - Zi]) 



X 
X 
X 
X 



xxx 
x x 
x 




[m 
[zi - z ,z 2 



x 

X 
X 
X 



Here the last column of the previous Q is discarded so that the remaining m = n 
columns form an orthonormal basis for the Z-flat tangent space, 7{Z). 



Z-flat gradient of expected return (rtndecomp: 204 223) 
g T F(:, 2 : n) = E(2 : n) — e\\ # eflag = false if solution is exact. 

e = e 1 - g T F(:, 1); e F = ||g||; # when eflag = false. 

If expected return is an affine function of (vector) risk, then 

[e 2 - ei, . . . , e n - ei] = g T [z 2 - z 1 , . . . z n - z x ] 

can be solved exactly for g G 7(Z). This will be the case (and eflag will be false) 
unless 1m is parallel to the i?-flat (Proposition 1). 

In this discussion we will assume that 1m is not parallel to the i?-flat. Then Ui = 
g/||g|| G 7(Z) is direction of steepest increase in expected return, and 

e = eo + ei?Ui T z (A. 2) 

holds for all z in the Z-flat, where = e\ — g T Zi and cf = ||g||- 

In this block we solve for the coordinates of g with respect to the current orthonormal 
basis, Q = [q 1; . . . ,q m ], for the tangent space 7{Z\ Then Cq and e F are computed 
from the coordinate representation. In the following block, g G R m denotes the 
coordinate vector corresponding to g G 7(Z). 



Householder reflection (rtndecomp: 225 235) 
H = I m -Pvv T ; /3 = 2/||v|| 2 ; ^g = <5 1 ||g||; 
F := HF; 

(still F ~ [zi - z , z 2 - zi, . . . , z n - zi]) 

F(:,j):=F(:,j)+F(:,l)- (j 1 //) 

(now F ~ [zi - z , z 2 - z , . . . , z„ - z ]) 

The first basis vector, q 1; of the current basis Q has coordinate vector Si, the first 
column of the m x m identity matrix I m . If H is the Househoulder reflection of R m 
that maps g to <5i||g||, then the first column of QH is the direction of maximum 



Householder 
S reflection H 
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increase in expected return in 7(Z). Thus, after the replacements Q := QH and 
F := HF, we still have F ~ [zi — zo, Z2 — Zj, . . . , z n — zi], but now times the first 
coordinates of F produces the corresponding changes in expected return: 

e F F(l, :) = [ei - e , e 2 - ei, . . . , e n - ei]. 

Then the replacements F(:,j) := F(:,j) + F(:, 1) (j = 2, . . . , n) result in 

F ~ [zi - z , z 2 - z , . . . , z n - z ] 

and 

£ = [ei, . . . , e n ] = e + e F F(l, :) 

with V = nf^ + F T F. 

Principal components of nonproductive risk (rtndecom: 237 243) 
F(2 : m, :) = [7S\^ T ; 7^ compact singular value decomposition. 

£ = diag(r2, . . . , r m ); # the 7* of Definition 1, Section 4.3. 

F(2 : m,:) := ZV T ; 
(final F ~ [zi - z , z 2 - z , . . . , z n - z ]) 

Entering this block the first row of F consists of the coordinates of [z\ — zq, z 2 — 
zo, . . . , z n — zo] in the productive risk direction. The remaining m — 1 rows represent 
the nonproductive risks. 

The principal components of nonproductive risk are computed with the code 

[U, S, V] = svd(F(2 : m, :), 0); 
The output matrices U and V (not to be confused with the covariance matrix V) have 
orthonormal columns and dimensions (m — 1) x (m — 1) and n x (m — 1), respectively. 
The principal nonproductive risks, r 2 > . . . > r n > 0, of Definition 1 are contained in 
the diagonal matrix S ( = S in the rtndecomp code). Now the replacements 

F{2 : m, :) := S V T , Q(:, 2 : m) := Q(:, 2:m)U 
maintain F ~ [zj — z , z 2 — z , . . . , z n — z ], but organize the last m — 1 rows of F 
in the principal directions of nonproductive risk, with row norms \\F(i, :)|| = r, for 
i — 2, . . . ,m. 

In this discussion we have described how the basis matrix Q changes from one block 
to another. It is always the case that Q T Q = I, though the size of the identity matrix 
diminishes from nxiias the algorithm progresses. The actual rtndecomp algorithm 
makes no attempt to keep track of the changing Q. 
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B The Octave algorithm 

# Function: [E, F, fO , eO , eF ] = rtndecomp ( rtns , wgts , pput) 

# 

# Purpose 

# To decompose financial return data into orthogonal 

# "systemic", "productive", and "nonproductive" 

# risk— factor s . 

# Input 

rtns — M by n matrix of periodic returns . 

# wgts — M by 1 vector of positive weights or a scalar , 

# default: wgts = ones (M, 1) / M if wgts is a 

# scalar or rtns is the only input. 

# pput — periods per unit time, default: pput = 1. 



# example: pput = 252 market— day s / year . 

# Output 

# E — 1 by n matrix of expected returns . 

# F — m by n matrix of risk coefficients. 

# rank(F) — m unless F = zeros (1 , n). 

# fO — systemic risk ( nonnegative ) . 

# eO — systemic expected return. 

# eF — return per unit change in the F(l, :)— direction . 

# (eF >= 0) 

# Global output 

# eflag — true iff the constant return vector, ones (M, 1 ) , 

# parallels the returns flat or, said another way, 

# ones (M, 1 ) = rtns * x 

# for some n— vector x with sum(x) = 0. 

# Variables and relationships 

# The 1 x n expected return matrix is 

# E — ( wgts * pput ) ' * rtns . 

# The n x n covariance of return matrix is 

# V — Z' * diag(wgts * pput) * Z, 

# where Z is the M x n matrix of risk vectors 

# Z — rtns — ones (M, 1) * E / pput. 

# The output variables satisfy 

# 1) V = f0"2 + F' * F. 

# 2) E — eO + eF * F(l, :) unless eflag is true; 

# then this relationship is approximate . However 

# mean(E) = meanfeO + eF * F(l, :)) is always true. 

# We refer to norm(F(l, :)) as the productive risk 

# when eF is nonzero . 

# 3) The row norms, 

# tau (i) = norm(F( i , : ) ) ( i — iO , . . , m) , 
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# are the principal nonproductive risks , where 

# iO = 2 if eF > 0, and iO = 1 if eF = 0. 

# The nonproductive risks tau ( i ) decrease with 

# increasing i, and the corresponding row vectors 

# are pairwise orthogonal: 

# F(i, :) ' * F(j, :) = (iO <= i < j <= m) . 

# 

# 

# Copyright (C) 2012 Vic Norton <mailto:vic@norton.name> 

# 

ft This GNU Octave program 

# http :// www. gnu . org / s oftware / octave/ 

# is free software : you can redistribute it and/or modify it 

# under the terms of the GNU General Public License 

# http : / / www. gnu . org / copyleft / gpl . html 

# as published by the Free S oftware Foundation 

# http://www.fsf.org/ 



# This program is distributed in the hope that it will be 

# useful, but WITHOUT ANY WARRANTY; without even the implied 

# warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 

# PURPOSE. See the GNU General Public License for more 

# details . 



function [E, F, fO , eO , eF] = rtndecomp ( rtns , wgts , pput ) 

if !( nargin >= 1 nargin <= 3) 

usage(" [E, W F , w f0 , w e0 , u eF] ^= w rtndecomp (rtns , ^ wgts , w pput ) " ) ; 
endif 

[M, n] = size(rtns); 
global eflag ; 

epsO — eps * le2 ; # precision 

## check/set wgts & pput 

if nargin = 1 || i s s c al ar ( wgts ) 

wgts = ones(M, 1) / M; 
else 

if ! is vect or ( wgts ) 

error ( " wgts^must^be^a^ scalar „or„a„vector" ) ; 
endif 

if length(wgts) != M 

error ( " wgts ^must w have ^length JYI" ) ; 
endif 
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ss if any(wgts <= 0) 

89 error ( " wgts w must^be^ positive " ) ; 

90 endif 

91 if rows(wgts) = 1 

92 wgts = wgts ' ; 

93 endif 

94 wgts /= sum(wgts); 

95 endif 

96 if nargin = 3 

97 if pput < 1 

98 error ( " g^must^be^ 1 ^ or ^greater 11 ) ; 

99 endif 

100 else 

101 pput = 1 ; 

102 endif 

103 sqrtpput = sqrt(pput); 

104 

105 /t4^ expected return matrix E and "risk" vectors Z for pput — 1 

106 E = wgts ' * rtns ; 

107 if nargout <= 1; 

108 E *= pput ; 

109 return; 
no endif 

in sqrtwgts = sqrt(wgts); 

112 Z = diag ( sqrtwgts ) * rtns; 

n3 Z — = sqrtwgts *E; ft V = Z' * Z = covariance matrix 

n4 sigZ = norm(Z, "fro"); ft square root of total variance 

115 

n6 ## QR factorization of Z with column pivoting 

n7 [Q, F, J] = qr(Z, 0); # Q is not used 

us m = rows (F ) ; 

no epsZ = epsO * sigZ ; 

120 while m > 1 && norm(F(m, m : n)) <= epsZ 
F(m, :) = []; 

122 m — = 1 ; 

123 endwhile 

124 Jinv = [1 : n] * eye(n)(: , J) '; 

12.", 

126 ## finish up when m = 1 and all co efficients of F are the same 

127 if n = 1 

128 fO = abs(F) ; 

129 if pput > 1 

130 E *= pput ; fO *= sqrtpput ; 

131 endif 

132 eO = E; eF = 0; F = 0; eflag = 0; 
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133 return; 

134 else if m = 1 

135 Fmin = min(F); Fmax = max(F); 

136 if Fmax — Fmin < epsO * max (abs( Fmin ) , abs(Fmax)) 

137 eO = mean(E ) ; 

138 fO = abs((Fmin + Fmax)/2); 

139 eF = ; 

140 F = zeros ( 1 , n ) ; 

141 ee = sum(E . * E); 

142 eflag = ee — n * eO * eO > epsO * ee ; 

143 if pput > 1 

144 E *= pput ; eO *= pput ; fO *= sqrtpput ; 

145 endif 

146 return; 

147 endif 

148 endif 

149 

150 /t4^ A-t this point the columns of F are the coordinate vectors 

151 # of the columns of Y — [yl , y2 , ■ ■ ■ , yn\ — Z(: , J) with 

152 # respect an orthonormal basis for the linear space, L(Z) , 

153 # spanned by the columns of Y or Z. We represent this 

154 # situation with the notation 

155 # F ~ [yl , y2 , . . . , yn] . 

156 # As the algorithm progesses the orthonomal basis for L(Z) 

157 # changes. We do not keep track of the changing basis , but 

158 # we do continue to note the elements of L(Z) represented by 

159 # the coordinate — vector columns of F. 

160 

161 ## E—flat & Z—flat tangent spaces 

162 el = E(J(1)); # base expected return 

163 yl = Z(: , J ( 1 ) ) ; # base risk vector 

164 nml = n — 1 ; 

165 B = E(J(2 : n)) — el; # differential expected returns 

166 F(l, 2 : n) — — F(l, 1); # differential risks 

167 

168 77^ QR decomposition of Z—flat tangent space by Givens 

169 # rotations 

170 for j = 2 : m 

171 jml = j - 1; 

172 [cs, sn] = givens (F(jml , j), F(j , j)); 

173 GV = [ cs , sn ; — sn , cs ] ; # Givens rotation 
F(jml : j , 1) = GV * F(jml : j , 1); 

175 F(jml : j , j : n) = GV * F(jml : j , j : n); 

F(j , j) - 0; 

177 endfor 
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178 # F [yl , y2 - yl , . . . , yn - yl] 

179 

180 ## systemic risk fO 

181 if norm(F(m, m + 1 : n)) <= epsZ 

182 ## the Z— flat tangent space has dimension m — 1 
iss fO = abs(F(m, 1)); 

F(m, :) []; 
iss m — = 1 ; 

186 else 

187 the Z— flat tangent space has dimension m 
iss fO = 0; 

189 endif 

wo # F ~ fyl - yO , y2 - yl , . . . , yn - yl ] 

191 # where yO = zO is the point in the Z—flat that is closest 

192 # to the origin . 

193 

194 check for constant E 

195 Emin = min(E); Emax = max(E); 

196 if Emax — Emin < epsO * max(abs(Emin) , abs(Emax)) 

197 E is constant : finish up 

198 F(: , 2 : n) += F(: , 1) * ones(l , n - 1) ; 

199 # F ~ fyl — yO , y2 — yO , ... , yn - yO] 

200 [U, S, V] = svd(F(: , Jinv) , 0); 

201 F = S * V ; 

202 E *= pput ; eO = mean(E); 

203 eF = 0; fO *= sqrtpput ; F *= sqrtpput ; 

204 return; 

205 endif 

206 

207 ## Z—flat direction of maximum increase 

208 # in expected return 

209 g = (B / F(: , 2 : n)) '; # E— gradient if eflag 

210 eflag = norm(g' * F(: , 2 : n) — B)/norm(B) > epsO ; 

211 if ! eflag # E = eO + eF * F(l, :) will be exact 

212 eO = el - g' * F(: , 1); 

213 eF = norm(g ) ; 

214 # F ~ fyl — yO , y2 — yl , ... , yn — yl] 

215 else # E — eO + eF * F(l, :) will be approximate 

216 F(: , 2 : n) += F(: , 1) * ones(l , n - 1) ; 

217 # F ~ fyl — yO , y2 — yO , ... , yn — yO] 

218 v = mean(F') '; # v — ymean — yO 

219 G = F — v * ones ( 1 , n ) ; 

220 7^ G ~ fyl — ymean, y2 — ymean, ... , yn — ymean] 

221 emean = mean(E); 

222 C = E(J) — emean; 
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223 g = (C / G) '; # E— gradient if ! eflag 

224 eO = emean — g' * v; 

225 eF = norm ( g ) ; 

226 endif 

227 

228 the Householder reflection , H g — (x, 0, ... , 0) 

229 # puts the productive risk coordinates in F(l, :) 

230 [v, b] = housh(g, 1, 0); 

231 signx = sign((g — (b * v) * (v' * g))(l)); 

232 F=F— (b*v)*(v'*F); 

233 F ( 1 , : ) *= signx ; 

234 if ~ eflag 

235 F ~ [yl — yO , y2 — yl , ... , yn — yl] 

236 F(: , 2 : n) += F(: , 1) * ones(l , n - 1) ; 

237 endif 

238 # F ~ fyl — yO , y2 — yO , ... , yn — yO] 

239 

240 P u t th e principal components of nonproductive risk 

2-11 # in F(2 : m, :) 

242 if m > 1 

243 [U, S, V] - svd(F(2 : m, :) , 0); 

244 F(2 : m, :) = S * V ; 

245 endif 

246 # F ~ fyl - yO , y2 — yO , ... , yn - yO] 

247 

248 unpermute columns of F 

249 F = F ( : , Jinv ) ; 

250 # F fzl - zO , z2 - zO , ... , zn - zO] 

251 

252 ## scale these results if pput !— 1 

253 if pput != 1 

254 E *= pput ; eO *= pput ; 

255 eF *= sqrtpput ; fO *= sqrtpput ; F *= sqrtpput ; 

256 endif 

257 

258 endfunction 
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